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ABSTRACT 

This thesis studies the estimation from interarrivai and 
service time data of the mean number of customers in service 
at time t for an M/G/6o queue. Two situations are considered. 
In one the parametric form of the service time distribution 
is known. In the special case in which the service time 
distribution is exponential the approximate bias and vari- 
ance of the estimate are derived and simulation is used to 
study an approximate normal confidence interval procedure. 
Simulation is also used to illustrate that assuming a wrong 
parametric model can lead to misleading results. In the 
other situation, the parametric form of the service time 
distribution is unknown and the empirical distribution of 
the service times is used in the estimate of the mean number 
of customers in service. In the case in which the customer 
arrival rate is known the distribution of the estimate is 
derived and an approximate normal confidence interval proce- 
dure is suggested. The use of the bootstrap and jackknife 
procedure to estimate variability and construct confidence 
intervals for the estimate is also studied both analytically 
and by simulation. 
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I. INTRODUCTION 



A. DESCRIPTION OF THE PROBLEM 

The application of probability theory to a wide variety 
of congestion problems has been described in many papers and 
books [Refs. 1,2,3]. Results of queueing theory are 
presented in terms of component distribution functions and 
stochastic processes (renewal, Poisson, etc) that are taken 
as known; only rarely are issues addressed that arise when 
actual data is to be used as a basis for inference from the 
models; however, see Cox(1965) [Ref. 4]. 

The concern of this thesis is inference problems for a 
particularly simple queueing model, the M/G/oo queue. In 
this model, customers arrive according to a Poisson process 
with rate and there are an unlimited number of 

independent servers. Service times for each server are 
independent, identically distributed with distribution 
function F. Let X(t) be the number of customers being 
served at time t. It is well known that if there are no 
customers being served at time 0, then 

P[X(t)=n] = exp[-M(t)j (1.1) 

k j 

where 

M( t) = X (^'^F( s)ds 

with F(t)=l-F(t) [Ref. 2] . Thus the distribution of X(t) is 
Poisson and is characterized by its mean M(t). 

In this thesis we will assume that the service time 
distribution F(t) is unknown and must be estimated from 
service time data and that the arrival process is known to 
be Poisson, except possibly for its rate . We will study 
the estimation of the mean number of customers being served 
at time t, M(t). 
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B. SCOPE OF THE THESIS 



The purpose of this thesis is to study the estimate of 
the mean number of customers being served at time t for a 
M/G/oo queue. This mean completely characterizes the 
distribution of the number of customers being served at time 
t. We will assume that the service time distribution and 
possibly the customer arrival rate are unknown and must be 
estimated from data. 

We generally divide the estimation method into two cases 
which we shall call "parametric estimation" and 
"nonparametric estimation". In the parametric estimation 
case, a particular probabilistic model is specified for the 
service time distribution and the parameters of the 
distribution are estimated. The resulting estimate of the 
survivor function is then used in the estimate of the 
expected number of customers being served at time t. In the 
nonparametric estimation method, the empirical survivor 
function is used in the estimate of the expected number of 
customers. 

In most cases, parametric assumptions concerning the 
service time distribution are difficult to justify. Hence 
nonparametric estimation procedure may well be preferred to 
parametric estimation when actual data is used. However, 
the nonparametric estimates can be expected to be less 
efficient than the parametric ones. 

The thesis is organized as follows. In Chapter II, the 
transient distribution for the number of customers being 
served at time t for the M/G/CD model is described and the 
equilibrium distribution as time goes to infinity is 
obtained. In Chapter III, we study parametric estimates of 
the mean number of customers being served under several 
assumptions for service time distributions. In the special 
case in which the service time distribution is exponential 
the approximate bias and variance of the estimate are 
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derived and simulation is used to study an approximate 
normal confidence interval procedure. Parametric estimates 
for gamma, mixed exponential, and lognormal distributions 
are also considered. Simulation is used to study the effect 
of assuming a wrong parametric model. In Chapter IV, a 
nonparametric estimate of the mean number of customers being 
served is described. This estimate is based on the 
empirical distribution of the service times. In the case in 
which the customer arrival rate is assumed known the 
distribution of the nonparametric estimate is derived and an 
asymptotic normal confidence interval procedure is 
suggested. The jackknife and Bootstrap methods for 
obtaining approximate confidence intervals are also 
described. The different estimators are compared by 
simulation. Chapter V describes the simulation and gives 
the results. 

In summary, this thesis studies the use of estimates of 
the service time distribution to obtain estimates of the 
mean number of customers being served for a M/G/X) queue. 
Both parametric and nonparametric estimates are considered 
and compared by simulation. 
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II. M/G/ OO QUEUE MODEL 



The M/G/oo queueing model is specified by the following 
assumptions. There are infinite number of servers. 
Customers arrive for service according to a Poisson process 
with rate \ . Service times are nonnegative independent 
identically distributed random variables with distribution 
function F. When a customer arrives, he immediately starts 
service. 

Let X(t) represent the number of ustomers in service at 
time t. It is well known that if chere are no customers 
being served at time 0, then 



P{X( t)=ki 



”"k.T” 



exp[->\p( t) ] 



( 2 . 1 ) 



where p(t)= [ 1-F( s ) ] ds: that is, X( t) has a Poisson 

distribution with mean \p(t) [Ref. 2]. Taking the limit as 
t -^00 in equation 2.1, we obtain the equilibrium 
distribution 



lim P[X(t)=kj = exp[ - A 1-F( X )dx] (2.2) 

t-.cc k/ 



Thus, the limiting distribution of X(t) as t-^oo is also 

y 00 

Poisson with mean ^ m, where m= F(x)dx is the mean service 
time. Therefore, the distribution of the number of 
customers being served at time t is Poisson with mean 



M( t) = A [Jf( x)dx 



(2.3) 



Here, the distribution of the number ' of customers being 
served at time t is characterized by value of its mean M(t). 
The value of M( t) depends upon the service time distribution 
which is assumed unknown and must be estimated from data. 



12 



This thesis considers the problem of estimating M(t) 
service and interarrival time data. 
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III. PARAMETRIC ESTIMATION METHOD 



A. DESCRIPTION 

In this chapter, it will be assumed that the parametric 
form of the service time distribution is known. In this 
case the estimation of the mean number of customers being 
served at time t, M(t), can be considered to be a function 
of the parameter estimates of the distribution. In 
particular, the estimate of M(t), when a parametric form of 
the service time distribution is assumed, is denoted by 
Mp( t ) , then 

Mp(t) =>\ C F(s)ds (3. 1) 

where F( t) is a survivor distribution of an assumed 
parametric form. 

In this chapter the rate of the arrival process will be 
assumed to be unknown. Maximum likelihood estimates of the 
mean interarrival times and the mean service times are used 
in the estimate of Mp(t). 

Four parametric service time distributions will be 
considered: the exponential, the gamma, the mixed 
exponential, and the lognormal distribution. In the 
exponential case, moment approximations are used to assess 
the bias of the estimate and to develop a confidence 
interval procedure based on asymptotic normality. The 
performance of the confidence interval procedure is assessed 
by simulation. 

In the remaining three parametric models, simulation is 
used to assess the performance of the parametric estimates. 
Another source of error in using a parametric estimator is 
that the wrong parametric form may be used. The effect of 
using the (wrong) exponential model in these cases is also 
assessed by simulation. 



14 



Each simulation has 300 replications; each replication 
consists of 50 independent service times from the specified 
distribution, and 50 independent interarrival times from an 
exponential distribution. The average relative bias and the 
average relative square error of Mp(t) are used to evaluate 
the performance of the parametric estimation method. All 
simulations were carried out on an IBM 3033 computer at the 
Naval Postgraduate School using the LLRANDOMII, random 
number generating package [ Ref. 6] . 

B. EXPONENTIAL SERVICE TIME 

In this section it will be assumed that the service time 
distribution is exponential; that is, F(t) = 1-EXP( -t/^. ), 
where it is an unknown parameter and must be estimated from 
the observed data. The maximum likelihood estimate of is 

A . 'n lu 

^ X;, where x- is the service time of the i customer. 
We will also assume that the rate of the Poisson arrival 
process \ is unknown and must be estimated. 

The interarrival times of the customers are denoted by y, 
/ Yi. / • • / y<*v • Since the arrival process is Poisson with rate 
\ , the interarrival times are mutually independent, 
positive random variables with the exponential distribution 
function having mean . The maximum likelihood estimate 
of \ is =n/ ^ y; . For an exponential service time 
distribution, an estimate of the mean number of customers in 
service at time t for an M/G/D queue is 

Mp(t) =\-M ( l-exp[ -t/u] ) (3.2) 

The estimate is a nonlinear function of the estimated 

A A 

parameters, \ and -U . In most cases, when estimating a 
function of the estimated parameters, bias is created by the 
nonlinear relationship of the estimated parameters. 

A 

Approximate formulas for the bias and variance of Mp(t) will 
now be derived. 
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\ A 1 

Let 3 be the mean service time, (3 = — £ , i=l , 2 , . . , n, 

and cs( be the mean interarrival time, ^ = — s y. , i=l,2,..,n. 

on Vi A 

By assumption, and Y;, are independent. The estimate of 

A 

Mp(t) can be represented by a function of the parameters o< 

A 

and 8 as fbllows: 






(3.3) 



There are no simple, exact formulas for the mean and 
variance of the quotient of two random variables. However, 
there are approximate formulas which are sometimes useful. 
The approximation can be obtained from the partial Taylor 
series expansions of M(o(,^) about the true means, o(. and (3 . 
The expansion is 



I ^ A 



(3.4) 



Since we assumed that the arrival process and the service 
times are independent, the covariance terms turn out to be 
zero when we take the expectation of both sides of equation 
3.4. Thus, we get ^ 

E[M(S^,^)] = M(o/,| 3) + ^ -'^-M(o{,3)Var(<i) + ^ ) Var(^ ) 



+ Ri 



(3.5) 



where converges to zero at the rate — ^ . The variance of 



estimate is 

-> A 



Var[M(o<,^)] = [ -^-M( c/ ,(3 ) J ^Var(3() + [ -^-M(0^ ,^ ) ] 2Var(^ ) 



+ R 



n 



(3. 6) 
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with R tending to zero at the rate . An approximate 
bias term, denoted by (9p(t), can be derived immediately from 
the equation 3.5, that is , (3p( t ) =E[ M(o^ , ^ ) -E[ M( / ^ ) ] 1 • 

Subtracting ^p(t) from the parametric value to correct the 
bias, leads to the bias corrected estimate of Mp(t). 

In order to compare the two estimates, bias and 
bias-corrected, we define the following notation. Let be 
the fraction of bias of Mp(t) against the true value M(t), 
and 0i be the fraction of square error of Mp(t) against the 
square of the true value: that is, 

e,= wh:.f3£« (3.7) 

( 3 . 8 ) 

where M( t) = \^^F( s)ds. 




Figure 3.1 Average Relative Bias of Mp(t) 
for Exponential with M =2 at t=l 



17 



A simulation experiment was performed to assess the 
performance of the estimates. In the i ^ replication, 50 
exponential interarrival times having mean 1 and 50 
exponential service times having mean 2 were simulated and 
estimates 



r^p(t) =\ M ( l-exp[ -t/u] ) (3.9) 

and 

Mp( t) = Mp(t) - t) (3. 10) 



where 

(3p(t) = ^ 5|rM(o(,^)Var(5^) + ^ ,^ ) Var(,S ) 

were computed. The estimated values of ^ and ^ were used 
in the variance formulas. The simulation was replicated 300 
times and the average relative biases 



e,(t) 



300 ; = i iMCt) 



(3.11) 



t) 



300 



I 



300 



A cl 



(3. 12) 



and the average square error 



6^( t) 



e"(t) 



_ j__ 


300 

z. 


j hd)- Mpd) j p 


( 3. 13) 


300 




Mlt3 






300 

z. 




(3. 14) 


300 


A=l 







were computed. 

Results of the simulation are presented in Figures 3. 1 
and 3.2. In Figure 3.1, the dotted line shows the average 
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Figure 3.2 Average Relative 3quare Error of Mp(t) 
for Exponential with U. =2 at t=l. ( n=50 , r=300 ) 

relative bias, 0^(t), as a function of t for the original 

A _c 

estimate Mp(t). The solid line shows 6,( t) for the bias 
^ c 

corrected estimate Mp(t). This superimposed figure 

indicates that 9\(t) for the bias-corrected estimate (with 

solid line) is almost constant and is small. The bias 

estimate produces large negative value of 9v( t) but ©\(t) 

approaches a limiting value as t yoD . Figure 3.2 shows of 

— ~ c 

the average relative square error 0j,( t) and 0i( t) plotted as 
a function of time. The dotted line gives t) and the 

solid line is 6^(t). It appears from figures that the 
estimate of bias described in equation 3.5 does correct for 
the bias. However, in Figure 3. 2 the bias-corrected 
estimate has a slightly higher relative square error than 
the original estimate. This higher relative square error 
could be due to correlation between the estimate itself and 
the estimate of its bias. 

Simulation was used to assess the performance of the 
following confidence interval procedure. In the i^*^ 
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TABLE I 



COVERAGE AND LENGTH OF 100( l-o()% C. I. FOR 
THE ORIGINAL ESTIMATE ( N = 50, R = 300 ) 



trial 


68 ? 




80 ? 




90 ? 




Length 

(s.d) 


C. R. 


Length 

(s.d) 


C. R. 


Length 

(s.d) 


C. R. 


1 


0. 2234 


10. 67 
72. 00 


0. 2879 


4. 00 
83. 00 


0. 3695 


1. 00 
90. 00 




( 0. 0415) 


17. 33 


( 0. 0535) 


13. 00 


( 0. 0687) 


9. 00 


2 


0. 2237 


9. 33 
73. 33 


0. 2882 


4. 33 
81. 00 


0. 3699 


1. 33 
86. 33 




( 0. 0396) 


17. 33 


( 0. 0510) 


14. 67 


( 0. 0655) 


12. 33 


3 


0. 2224 


11. 00 
68. 67 


0. 2866 


4. 33 
79. 00 


0. 3678 


1. 33 
87. 33 




( 0. 0410) 


20. 33 


( 0. 0529) 


16. 67 


( 0. 0678) 


11. 33 


4 


0. 2212 


9. 00 
66. 67 


0. 2851 


3. 67 
76. 00 


0. 3659 


1. 00 
83. 33 




( 0. 0402 ) 


24. 33 


( 0. 0519) 


20. 33 


( 0. 0666) 


15. 67 


5 


0. 2207 


13. 33 
62. 33 


0. 2844 


4. 67 
77. 67 


0. 3651 


0. 67 
88. 67 




( 0. 0409 ) 


24. 33 


( 0. 0527) 


17. 67 


( 0. 0677) 


10. 67 


6 


0. 2269 


11. 33 
70. 00 


0. 2925 


3. 67 
82. 67 


0. 3754 


0. 33 
89. 33 




( 0. 0433 ) 


18. 67 


( 0. 0558) 


13. 67 


( 0. 0717) 


10. 33 


7 


0. 2223 


10. 67 
69. 00 


0. 2865 


3. 00 
82. 00 


0. 3677 


0. 33 
89. 67 




( 0. 0375) 


20. 33 


( 0. 0483 ) 


15. 00 


( 0. 0620) 


10. 00 


8 


0. 2246 


7. 67 
71. 67 


0. 2895 


2. 33 
81. 33 


0. 3715 


0. 33 
87. 67 




{ 0. 0425) 


20. 67 


( 0. 0548) 


16. 33 


( 0. 0704) 


12. 00 


9 


0. 2232 


11. 00 
63. 67 


0. 2876 


6. 00 
82. 67 


0. 3692 


0. 67 
83. 67 




( 0. 0423 ) 


25. 33 


( 0. 0545) 


21. 33 


( 0. 0700) 


15. 67 


10 


0. 2270 


13. 67 
67. 67 


0. 2926 


3. 00 
81. 00 


0. 3755 


0. 33 
87. 33 




( 0. 0411) 


18. 67 


( 0. 0530) 


16. 00 


( 0. 0680) 


12. 33 


Average 


0. 2235 


10. 77 
68. 50 


0. 2881 


3. 90 
80. 63 


0. 3697 


0. 74 
87. 33 


( 0. 0410) 


20. 73 


( 0. 0528) 


15. 47 


( 0. 0678) 


11. 93 



replication of the simulation, 50 exponential interarrival 
times having mean 1, and 50 exponential service times having 
mean 2 were generated, and Mp( t) and Mp(t) were computed. 
The approximate variance of Mp(t) was computed for t=l using 
equation 3.6 with Rr\=0. The 100(l-o()% confidence limits L 
and U were computed by 
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TABLE II 



COVERAGE AND LENGTH OF 100( 1-(X)% C. I. FOR 
THE BIAS-CORRECTED ESTIMATE ( N=50, R=300 ) 



trial 


68 




80 ^ 




90 ? 




Length 

(s.d) 


C. R. 


Length 

(s.d) 


C. R. 


Length 
( s. d) 


C. R. 


1 


0. 2234 


14. 67 
69. 00 


0. 2879 


6. 67 
80. 67 


0. 3695 


2. 33 
89. 00 




( 0. 0415) 


16. 33 


( 0. 0535) 


12. 67 


( 0. 0687) 


8. 67 


2 


0. 2237 


13. 67 
69. 67 


0. 2882 


5. 67 
80. 33 


0. 3699 


2. 67 
86. 33 




( 0. 0396) 


16. 67 


( 0. 0510) 


14. 00 


( 0. 0655) 


11. 00 


3 


0. 2224 


16. 00 

66. 00 


0. 2866 


7. 67 
76. 33 


0. 3678 


2. 00 
87. 33 




( 0. 0410) 


18. 00 


(0. 0529) 


16. 00 


( 0. 0678) 


10. 67 


4 


0. 2212 


14. 33 
62. 33 


0. 2851 


6. 00 
74. 33 


0. 3659 


1. 67 
83. 00 




( 0. 0402) 


23. 33 


( 0. 0519) 


19. 67 


( 0. 0666) 


15. 33 


5 


0. 2207 


17. 67 
59. 67 


0. 2844 


9. 00 
74. 67 


0. 3651 


1. 67 
88. 67 




( 0. 0409) 


22. 67 


( 0. 0527) 


16. 33 


( 0. 0676) 


9. 67 


6 


0. 2269 


18. 67 
64. 00 


0. 2925 


6. 33 
81. 33 


0. 3754 


2. 00 

88. 00 




( 0. 0433) 


17. 33 


( 0. 0558) 


12. 33 


( 0. 0717) 


10. 00 


7 


0. 2223 


14. 67 
67. 00 


0. 2865 


5. 67 
80. 33 


0. 3677 


2. 00 
88. 67 




( 0. 0375) 


18. 33 


( 0. 0483) 


14. 00 


( 0. 0620) 


9. 33 


8 


0. 2246 


14. 00 
67. 33 


0. 2895 


4. 67 
80. 00 


0. 3715 


0. 67 
88. 00 




( 0. 0425) 


18. 67 


( 0. 0548) 


15. 33 


( 0. 0704) 


11. 33 


9 


0. 2232 


15. 67 
60. 00 


0. 2876 


9. 00 
70. 67 


0. 3692 


1. 67 
83. 33 




( 0. 0423) 


24. 33 


( 0. 0545) 


20. 33 


( 0. 0670) 


15. 00 


10 


0. 2270 


18. 33 
63. 00 


0. 2926 


9. 33 
75. 33 


0. 3755 


1. 67 
86. 00 




( 0. 0411) 


18. 67 


( 0. 0530) 


15. 33 


( 0. 0680) 


12. 33 


Average 


0. 2235 


15. 77 
64. 80 


0. 2881 


7. 00 
77. 40 


0. 3697 


1. 84 
86. 83 


( 0. 0410) 


19. 43 


( 0. 0528) 


15. 60 


( 0. 0678) 


11. 33 



L - Mp(t) - 



and 



U 



— Mp( t) + J 



(3. 15) 



(3. 16) 
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where z is the upper 1- ^ point of the standard normal 
distribution. Tables I and II show the results of 10 
independent simulations for the original and the 
bias-corrected estimate. Each simulation was replicated 300 
times. Tables report the average and standard deviation of 
the normal confidence interval length; the proportion, p of 
the intervals that covers the true value; the proportion of 
intervals that are too high, ( e. g. M( t) < L ); and the 
proportion of intervals that are too low,(e.g. M(t) > U ). 

Since the simulation replications are independent, it is 
possible to assess the uncertainty of p. If the confidence 
interval procedure is correct, then p should be within 



approximately ± 



of l-o(. 



The coverage rate in the 



3 oo 

tables indicate that the parametric estimates tend to 

A 

underestimated. Obviously, the distribution of Mp(t) is 
skewed right. However, the confidence interval procedure 

works well, regardless for both the original and the 

bias-corrected estimate. Both have a variable coverage 
rate. The difference of performance between two estimates 
is not significant. 



C. OTHER SERVICE TIME DISTRIBUTIONS 

1. Mixed Exponential Service Time 

In this subsection, service times having a mixed 
exponential parametric form will be considered. Customers 
arrive according to a Poisson process with unknown rate 
\ which must be estimated. Customers are of two types; with 
probability Pj , a customer's service time is exponential 
with mean ,u, ; with probability P^, , the customer's service 
time is exponential with mean M».. If T is the service time 
of an arbitrary customer, then 

P(T > t) = P, exp[ -t/M, ] + Pxexp[ -t/Ma.] (3.17) 
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where + P = 1. In this case, the mean number of 

customers being served at time t is 

Mp( t) = \ j P, M, Cl - expC-t/uO) c i - exp C-t/uO')j( 3 . 18) 

We will assume that a customer's type and service time are 
observable. 




Figure 3.3 Average Relative Bias of Mp(t) for 
Mixed Exponential with w,=2, xAu=. 75, P, =. 2 at t=l 

( n=50, r=300 ) 

A simulation experiment was done to assess the 
performance of Mp(t). In the i'^'^ replication, service times 
for 50 customers were generated, where the proportion P, of 
them have a type 1 and the proportion Pj_ have a type 2. A 
random number was drawn to determine the type of customer. 
If a customer was of type i, a service time was generated 
from the exponential distribution with mean ttc . For the 
simulation, the proportion P, is assumed known. The 50 
independent interarrival times were also generated. In this 
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Figure 3.4 Average Relative Square Error of Mp(t) for 
Mixed Exponential with u, =2, Ui.= . 75, Pj =. 2 at t=l 

( n=50, r=300 ) 

case , P(=0.2, =2. , M:>.=0.75, and \=1. The estimate Mp(t) 

was computed; the estimate of P^ is the proportion of 
customers that were type i; the estimate of the mean service 
time Ua. was the mean service time of type i customers; the 
estimate of \ was the mean interarrival time. The estimate 

A 

M 0 (t) of M(t) assuming that the service time distribution is 
exponential was also computed, that is, 

M^(t) = \- M ( l-exp[ -t/w] ) (3.19) 

with W equal to the mean service time and \ is equal to 
the mean interarrival time. The procedure was replicated 
300 times. Results of the simulation appear in Figures 3. 3 
and 3.4. 0»(t) is the average relative bias of the estimate 

and ^(t) is the average relative square error as computed 
in equation 3. 11 and 3. 13. The solid line is for the 

A 

correct parametric model estimate Mp(t). The dotted line 

A 

represents the exponential model estimate Mg(t). 
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In both figures, we compare the correct model 
(represented by a solid line) with the erroneous exponential 
model. As expected, the exponential model clearly does not 
perform as well as the mixed exponential model. In Figure 
3 . 3 , the level of bias for the exponential model is very 
high early in time, but is reduced, and stabilizes as t -^oo . 
Notice that the estimate using the exponential model appears 
too large. The average relative square error of the 
estimates is shown in Figure 3.4. Although the results of 
exponential model shows a slightly higher mean square error 
than that of the mixed exponential model, the results of 
exponential model are not too much worse. This is not 
surprising, since as t 4 C 0 the estimate just depends on the 
mean. 

2. Gamma Service Time 

In this subsection, the parametric form of the 
service time distribution is gamma. The probability density 
function is 



f(t) = -A- (3.20) 

where k and 3 are strictly positive parameters of the 
distribution, and k is further assumed to be an integer. By 
successive integrations by parts, we get 

F(t) = 1 - Z (3.21) 

A=o A y 



its mean and standard deviation are 



yU - 







JR 



Thus, k is the parameter that specifies the degree of 
variability of the service times relative to the mean. 

Since the arrival process is Poisson with rate \ , 
the interarrival times are mutually independent, positive 
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random variables with the distribution function 
G( y )=1-EXP( - A, y ) , where \ must be estimated from 
interarrival data y^ , i=l to n. Thus, the Mp(t) in this 
case is obtained by the successive integrations by parts of 
the survival function of the gamma distribution, that is 



r > ^ ;=o K I ^ 



( 3 . 22 ) 




Figure 3.5 Average Relative Bias of Mp(t) 
for Gamma with 3=1/ k=2 at t=l. (n=50, r=300 ) 

The performance of Mp (t) was evaluated by the 
simulation. In the simulation, the parameter k of the gamma 
distribution was assumed to be known, but the rate of the 
arrival process is unknown and is estimated. Two simulation 
cases were run. In the i'^'^ replication of the simulation, 50 
independent service times were generated, where the first 
simulation case used the gamma distribution having 3=1 and 
k=2 and the second case used the gamma distribution having 
3 =0. 5 and k=4; 50 independent interarrival times having 
3=1 were also generated. For k=2, the estimate is 
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( 3. 23 ) 



M_(t) = l-exp[ -^t] - t-exp[-At]} 

and for k=4, the estimate is 



Mp(t) = \ {- 7 -( l-exp[ -at] - exp[--at](3t + jS t^ 



AT) A "n 

where A =n/21 y- and /3 =kn/x x^ were calculated. 

\at\ ^ V-V 

based on an erroneous exponential service 
parametric estimator was also calculated 




Al 



^ t* 
( 3. 2 



11 



An estimate 
time model 



Me(t) 



= \ M ( l-exp[ -t/u] ) 



( 3. 25) 



with M = — Z X; 




Figure 3.6 Average Relative Bias of Mp(t) 
for Gamma with Q=. 5, k=4 at t=l. ( n=50, r=300 ) 

The simulation was replicated 300 times. The 
average relative bias 6,(t) and the average relative square 
error 9o.(t) were calculated. These results appear in 
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Figure 3.7 Average Relative Square Error of M (t) 
for Gamma with 3 =1/ k=2 at t=l. ( n=50, r=300 ) 

Figures 3.5, 3.6, 3.7, and 3.8. The dotted lines in the 
figures show the results of the erroneous exponential model. 
The solid line represents the results of the gamma model. 
The figures clearly show that the performance of the 
exponential model is not as good as that of the gamma model 
initially. However, both models have the same equlibrium 
state for t > 15, approximately. Figure 3.5 shows the 
average relative bias of the estimate in the case of k=2 and 
Figure 3. 6 shows the same in the case of k=4. In both 
figures, the erroneous exponential model has a high level of 
bias and the bias of the gamma model is almost constant; 
however, both models have exactly the same limiting value as 
t -^00 . The average relative square error appears in the 
Figures 3.7 and 3.8. Figure 3.7 represents the average 

A A 

relative square error of Mp(t) and Mg(t). in case of k=2, 
and Figure 3. 8 represents the same in the case of k=4. In 
Figure 3.7, the exponential model has a poorer performance 
than the gamma model, but the difference is small. In 
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Figure 3.8, the exponential model has a large value of mean 
square error initially and the level of mean square error 
associated with the exponential model grows as the value of 
parameter k is increased. 




Figure 3.8 Average Relative Square Error of Mp(t) 
for Gamma with 3 =. 5, k=4 at t=l. ( n=50, r=300 ) 

3. Lognormal Service Time 

In this subsection, the service time distribution is 
assumed to be lognormal. Let X be a random variable, and 
let a new random variable Y be defined as Y=ln X. If Y has 
a normal distribution, then X is said to have a lognormal 
distribution. The density of a lognormal distribution is 
given by 

f(x) = --l---exp[ - ----(In X -£ )^] (3.26) 

X. 

where -C0< § <00 and <f >0. Set Y=ln X - § and by the 



29 



integration 



= 




1 




4an<r 


J-00 


~ ( 


X - ^ 


' Y ^ 


<T 


• ) is 


the st 



<2\p C' /V i cT'-J 






( 3. 27) 



where G^( • ) is the standard normal distribution function. 
If the service time has a lognormal distribution, then the 
estimated mean number of customers being served at time t, 

A 

Mp(t), is obtained by integration-by-parts 



M 



p(t) = K {t[l-F(t)l + j^^sF(s)dsl 



( 3. 28) 



Substituting equation 3.26 for f(t) and equation 3.27 for 
F(t) in the equation 3.28, we obtain 



A 



(t) = (3.29) 



where G^( ■ ) is the standard normal distribution function. 

The lognormal distribution is positively skewed and the 
level of skewness depends upon the value of mean and 

variance of the distribution. If the value of mean is 

decreased but the variance is increased, then the shape of 
distribution tends to be more skewed and it approaches the 
shape of the exponential model. 

The performance of the parametric estimate was 

assessed by simulation. In each replication of the 

simulation, 50 independent lognormal service times and 50 
independent exponential interarrival times were generated. 
The simulation generated two sets of the service times. One 

set of service times is from the lognormal distrbution with 

§=0.193 and cf^=l , and the other set is from the lognormal 
distribution with £=0.568 and <£^ =0.5. To estimate the mean 
and variance from the data, the logarithm of the service 
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Figure 3. 9 
Lognormal with 



Average Relative Bias of Mp( t) for 
^ =. 193, 0^=1 at t=l. ( n=50, r=300 ) 




Figure 3.10 Average Relative Bias of lyip(t) for 
Lognormal with §=.568, (f^=. 5 at t=l. ( n=50, r=300 ) 

time data, t;=ln x; , i=l to n was computed. The mean and 
variance are expected by 
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^ I fV ' ^ 

I = J = 

Thus the estimate is 



I ' ‘ ^ -s 

h-l 



k 



»(t) = \\ t-ri-Gry^- :>(;.( 



i,i 



' 

-l-<r 






( 3. 30) 



where \ is the estimate of the arrival rate. An estimate 
based on an erroneous exponential model 



M(.(t) =X-M(1 - exp[ -t/A] ) (3.31) 



A. I ^ 

was also computed, where m = — X x^. 

^ Xn ^ 




Figure 3.11 Average Relative Square Error of Mp(t) 
for Lognormal with §=.193, (J^ = l at t=l. ( n=50, r=300 ) 

The simulation was replicated 300 times. Figure 3. 8 
shows the tendency of 0v(t), the average relative bias of 
Mp(t), for the correct model (shown with a solid line) and 
the erroneous exponential model ( shown with a dotted line) 
for the simulation with §=0.193 and cT^=l in Figure 3.9, and 
with §=0.568 and <f =0. 5 in Figure 3.10. 
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Figure 3.12 Average Relative Square Error of Mp(t) 
for Lognormal with. ^=.568, 5 at t=l. ( n=50, r=300 ) 

In both figures, the average relative bias of the 
exponential model is also large initially. As expected, the 
average relative bias of the exponential model in Figure 
3. 10 is larger than the results in Figure 3. 9. As t -^co , 
the exponential model shows better performance and has a 
limiting value. Figures 3. 11 and 3. 12 show the tendency of 
0jt), the average relative square error, for the correct 
model ( shown with a solid) and the erroneous exponential 

model (shown with a dotted line). For Figure 3.11 the 

parameters, § =0. 193 and d'^=l, are used to generate data. 
And for Figure 3.12 the parameters, |=0.568 and <f =0.5, are 
used. The value of average relative square error of the 
exponential model in Figure 3. 12 is also higher than the 
results in Figure 3. 11. 

D. SUMMARY 

The general conclusions of this chapter are that the 
parametric estimation method is a highly efficient for 
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obtaining estimates of M( t) whenever the correct assumption 
for the model is given. The structure of the estimate of 
M( t) is clearly biased since it is a nonlinear function of 
the estimated parameters for the service time distribution 
and the customer arrival rate. However, for the service 
time distributions considered the indications are that the 
bias is small. Hence the parametric estimation method 
performs very well, whether or not the estimate is corrected 
for bias, when the correct parametric form is used. However 
the performance of the parametric estimation is very poor 
when the wrong parametric model is used. For instance, the 
erroneous exponential model often has a high level of bias 
and mean- squared error. Notice that the exponential model 
converges to the same limiting value as the correct model as 
t -^oo in all the cases considered. This is because as t •♦•oo 
all estimates use the mean service time to estimate the 
integral of the servivor functions. 
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IV. NONPARAMETRIC ESTIMATION 



A. DESCRIPTION 

Nonparametric methods are statistical techniques which 
are applicable regardless of the form of the distribution 
function that the measurement comes from. In this chapter, 
these techniques will, for the most part, be based on the 
order statistics. 

Let x,,x^ denote a random sample from a CDF F, 
and let s^^ , s^^^^ , . . . , s^^j denote that corresponding order 
statistics. Then the sample CDF is defined by 



I 

F^(t) = ( number of s^^.^ less than or equal to t) 

n ^ -t 3 

For fixed time t, F,^(t) is a statistic since it is a 

A 

function of the sample. In fact, for fixed time t, F„(t) 
has the same distribution as that of the sample mean of a 
Bernoulli random variable. We know by the central limit 

A 

theorem that F„(t) is a asymptotically normally distributed 
with mean F(t) and variance ( -^ ) F( t ) [ 1-F( t ) ] . 

Recall that (in chapter II) the mean number of customers 
being served at time t, M(t), is a function of the arrival 
rate \ and the survivor function of the service times. 
Hence a nonparametric estimator, denoted by M^(t), can be 
represented by the estimated values of \ and F(t). The 
estimated survivor function of service time is 



F^(t) 




0 






i f 0 ^ t < s, 






if s $ t < s . 

O'-! > 



if t ^ s 



for i=l , 2 , . . , n-1 



35 



Now, using the fact that M^(t)= \\^F(s)ds, 
nonparametric esimate 



Mw(t) 



= \ t 



ti 



'n 






V-' > 



if 0 $ t < s 

if s ^ t < s 

if t ^ s 



we obtain a 



(4. 1) 



where \ is the estimated arrival rate. Note :hat the 

nonparmetric estimate has a limiting value as t -t-oo , that 

is, lim M,j(t)=\m where m is the mean service time. 

•l-^oo 

In this chapter, we will consider two different 
situations. In one case, we will assume that the arrival 
rate is known but the distribution of service times is 
unknown and must be estimated. Based on this assumption, 

A 

Mg(t) is expressed simply in terms of the order statistics 
of the service times as follows 







A 




(4.2) 



when s <t4 s,,, . In the Appendix A, we derive the 

distribution of M^(t) in this case. Its mean and variance 
are 



E[M^(t)] =\^V(s)ds 



(4. 3) 



Var[M^(t)] = !^^sF(ds) - [ 5^sF(ds)]2 + t2F(t)F(t) 

- 2tF(t)[ S^^sF(ds)] j (4.4) 

A 

Thus, M^(t) is an unbiased estimate of M| 4 (t), Further, as 

A 

the sample size n is increased, M^(t) is asymptotically 
normal. Thus an approximate normal 100(l-c<)% confidence 

A 

interval for M^(t) is given by 
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(4.5) 



M^(t) ± z I Var[M ( t) j 

where z (,4 is the upper 1 - ^ point of standard normal 
distribution and Var[M^^(t)] is given in Appendix A. In this 
chapter, we will also study the jackknife and the bootstrap 

A 

procedures for obtaining confidence intervals for M^(t) in 
the case in which \ is known. In the second case, we will 
assume that the arrival rate >s is also unknown and must be 
estimated. Then, a nonparametric estimate M^(t) is the 
product of two estimates, 

^ A 

A ^ \ Ks cr\-K 

M^)(t) =\[ — Xs + (4.6) 

-ry 

where and K is the number of service times that 

are less than or equal to t. There are no exact functional 
forms for the mean and variance of ( t) in this case. 
However, the jackknife and the bootstrap methods can be used 
to obtain confidence intervals. This will be described 
below. 

B. JACKKNIFE ESTIMATION METHOD 

In this section, we will study the jackknife procedure 
for obtaining a confidence interval for M^ (t). The 
jackknife was first introduced by Quenouille (1949) for the 
purpose of reducing the estimate bias, and the procedure was 
later utilized by Tukey (1958), to develop a general method 
for obtaining approximate confidence intervals [Refs. 7,8]. 

The basic idea of the jackknife estimation method is to 
assess the effect of each of the groups into which the data 
have been divided, not by the results for that group alone, 
but rather through the effect upon the body of data that 
results from omitting that group. The two bases of the 
jackknife are that we make the desired calculation for all 
the data, and then, after dividing the data into groups, we 
make the calculations for each of the slightly reduced 
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bodies of data obtained by leaving out just one of the 
groups. A special case of Jackknife estimation is called 
the "complete jackknife estimation", where the number of 
subgroups is n (the size of sample); the subgroup is 

obtained by deleting the i'^'^ observation; thus the size of 
each subgroup is n-1 [Ref. 9]. Attention will be restricted 
to complete jackknife estimation in this study. 

A ^ 

Let Mf^.^(t) be the estimated mean number of customers 
being served at time t on the portion of the sample that 

I I ^ 

omits the i ^ sample. Let (t) be the corresponding 

estimator for the entire sample and define the i"^*^ 
Pseudo-value by 

M,(t) = nM^„(t) - (n-l)M,',(t) (4.7) 

A A jL 

The jackknife estimate (t) and an estimate of its 

variance are given by 

A . ^ A 

M (t) = -- X MUt) (4.8) 



s 



V. 

o 






[M;( t) 



M3(t)] 2 



(4.9) 



Tukey (1958) proposes that the n estimated pseudo values be 
treated as approximately independent and identically 
distributed random variables [Ref. 9]. Hence, the statistic 



v/fT C - KcuiC^) > 



has an approximate t-di stribution with n-1 
freedom, which leads to the approximate 
confidence interval 



(4. 10) 
degrees of 
100(l-o( )% 



M^(t) ± t,.^S, 



(4. 11) 
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where t is the upper 1- ^ critical point of the 

t-distribution with n-1 degrees of freedom. The confidence 
interval given by equation 4. 11 is a function of the 

estimated variance. In the remainder of this section, we 
will describe several methods of implementing the confidence 
interval procedure. We will also obtain an analytic 

expression for the jackknife estimate and its variance 

estimate for the case in which the arrival rate \ is known. 

1. Jackknife Estimate with Known Arrival Rate 

In this subsection, the arrival rate is asssumed 
known. In this case the closed form expression for the 
jackknife estimate and its variance estimate can be derived. 

The nonparametric estimate of the mean number of 

A 

customers being served at time t, M^t), can be expressed in 
terms of the service times as follows: 



M^(t) 



A 

Xs . + ----t] 
;=\ (■*> n 



(4. 12) 



where ^ the order statistics of independent and 

identically distributed random quantities from the unknown 

A 

probability distribution F, and the variable K is the number 
of ' s which are less than t. This equation shows 

A 

immediately that IVI^( t) is the linear function of the order 
statistics of service times. 

The jackknife estimate is based on sequentially 
deleting point and recomputing the estimator. Removing 
point from data set gives a different empirical 



O 

probability distribution 

/ . • , 

the estimate. 



with mass 






and a 



corresponding recomputed value of 






In the jackknife process, the 






)f 

pseudo 



value is 
M;(t) = 



\t 
\ s 



(A) 



if i > 



if i $ 



A 

K 

A 

K 



(4. 13) 
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for the fixed time t. Accordingly the pseudo values M-(t) 
have just K+1 different values. The jackknife estimate is 



M^(t) 






V) 



m- K 

n 



t] 



( 4. 14) 



This result is exactly the same as the original estimate. 

A 

This is because the estimate M 3 -( t) is unbiased. In Appendix 
B, the jackknife variance estimate is derived as follows: 

+ t^F^(t)F4t)j (4.15) 

A 

where F^(t) is the sample survivor function. Comparing 

/T\ 

equation 4.15 with equation 4.3 we see that (~^“~)[Var {My 
(t)j] = E[ Var [M^( t ) ] ] . Thus the jackknife variance estimate 
tends to be conservative in the sense that its expectation 
is greater than the true variance of M^(t). We will now 
describe two selected procedures to obtain confidence 
intervals for the jackknife estimate. Tukey suggested that 
the statistic in equati 'n 4. 10 has an approximate 
t-di stribution with n-1 degrees of freedom, which leads to 
the approximate two-sided 100( l-c<)% confidence interval 

M^(t) ± t Var[M3( t)] (4.16) 

for M^j(t), where t is the upper 1-^ critical point of 
the t-di stribution with n-1 degrees of freedom. However, 

A 

the n estimate pseudo values have just K+1 different values. 
Hence, another possible procedure is to adjust the degrees 
of freedom of the t-di stribution, that is, subtract one from 

A 

the number of different pseudo values (K+1), and use the 
result as the degrees of freedom. The length of confidence 
interval generated using the adjusted degrees of freedom (K) 
is slightly wider than that generated using the usual 






A 
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degrees of freedom (n-1) and the coverage rate should be 
increased. 

2. Jackknife Estimate with Unknown Arrival Rate 

In this subsection, it will be assumed that the rate 
of Poisson arrival process is unknown and must also be 
estimated. The maximum likelihood estimate is 
where y; is the interarrival time between i and (i-1) ^ 
customers. A nonparametric estimate of mean number of 
customers being served at time t is given by 






A 









(4. 17) 



where S(.;;,)'s are the order statistics of independent, 
identically distributed random quantities from the unknown 
probability distribution F. It is assumed that the s and 

A 

' s are independent. The variable K is the number of S ' s 
which are less than t. The data consist of two independent 
random samples, 

S , , , . . . , - F and Y, , Y^. , . . . , Y„ ^ Q 

F and Q being two possibly different distribution on the 
real line with Q, the exponential distribution with mean . 
From equation 4.17, the estimate M^i( t) is the product of two 
estimates. One is the function of y,c , \ = and the 

other is the function of s^ , H(s) = '^ ^ s^ + t. there are 
many possible ways to perform a two-sample jackknife 
procedure. We will call one method the "paired sample 
jackknife" procedure. Since the size of both samples is the 
same, we make the one set of observations by pairing 
respective observations, that is, ( s ,, y^ ),( s^ , y^ ,...,( s^ , y^ 

A 

). As with the one-sample jackknife, we estimate the M<^,^(t) 

^ A 

for all the data, and then, we estimate t) based on the 

I L 

remaining data obtained by leaving out just the i pair. 

-tU ^ 

Thus the i ^ pseudo value Mx(t) is 



M,(t) = nM^,((t) - (n-l)M;.^(t) 



(4. 18) 
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and the jackknife estimate Mj(t) and variance estimates are 
given by 

A 1 ^ 

Mj(t) = -^.^M,(t) 



Sj = 



I J1 A 

X[M-( t) - M^( t)] 2 
'vv c^-o ^ a' ' ' 

on these statistics, 



Based on these statistics, an approximate two-sided 
( 1“*^)% confidence interval is given by 



100 



M (t) ± 



A 

t I* oi S 



i-i 



4 



(4. 19) 



where t is the upper 1- ^ point of t-distribution with 

n-1 degrees of freedom. A second method is called the 
"separated sample jackknife" procedure. Since we assumed 
that the X;'s and ' s are independent, we can perform the 
jackknife procedure separately for each sample, and then, 
estimates which combine jackknife estimates and the 
jackknife variance estimate can be computed. 

A, A 

Let M 3 |(t) be the jackknife estimate of \ and Vy be 
the jackknife variance estimate for Let M 3 _j^(t) be the 

f "t — * ^ 

jackknife estimate of F(s)ds and be its jackknife 

variance estimate. Then the combined jackknife estimate of 
is 

M^^(t) = M3,(t)-M3^(t) 



( 4. 20) 



and the combined jackknife variance estimate is 



A A ^ 

s,. = 



A A 






(4.21) 



The approximate two-sided 100 ( l-o() % confidence interval is 
given by 



Mj^(t) ± t ,.a /TrT 



( 4. 22) 
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where t is the upper 1-^ point of t-distribution with n-1 
degrees of freedom. 

C. BOOTSTRAP ESTIMATION METHOD 

Efron(1979) introduced the bootstrap method for 
estimating the distribution of a statistic computed from 
observations [Ref. 10]. The bootstrap estimate is obtained 
by replacing the unknown distribution by the empirical 
distribution of the data in the definition of the 
statistical function. In practice, the distribution of the 
statistic is approximated by Monte Carlo methods. 

For convenience, the arrival rate is assumed to be known 
and equal to 1, then the nonparametric estimate M^(t) is 
just a function of service times. This is a one-sample 
problem. The bootstrap procedure is as follows: 

1. Suppose that the data points Xi ,Xi_ ,... ,Xr> are 
independent observations from the unknown distribution 
F. Then the true estimate is 

M^(t) = S^^F(s)ds (4.23) 



2. We can estimate the distribution F by the empirical 
probability distribution Fn . 

F : mass on each observed data point X; 
i=l,2, ,n 

3. The bootstrap estimate of M^j(t) is 



Ma(t) = Fn(s)ds 



( 4. 24) 



To obtain an estimate of variability for MQ(t), we precede 
as follows 

(1) Construct F^^ < the empirical distribution function, 
as just described. 

(2) Draw a bootstrap sample x^ , x^ , . . . , x^ by 
independent random sampling from F,, . 

Notice that we are not getting a permutation distribution 

since the values of are selected with replacement from 



43 



the set ( X, , , . . . , ). As a point of comparison, the 
ordinary jackknife can be thought of as drawing samples of 
size n-1 without replacement. 

(3) Compute an estimate of Mm (t) for each bootstrap 
replication, MwCt), that is, the value of statistic 
evaluated for the hootstrap sample. 



* -'-In-|l(x?<t)lt (4.25) 



where I(x^t)=\l if x^t 

\0 otherwise 

(4) Do step (2) some large number "B" times obtainin' 
independent bootstrap replications M'^' (t), 1 



Based on the bootstrap replications, the approximate 
estimate of M^j( t) and its variance are obtained by 






(t) 



(4.26) 



vir[M^( t)] 



I 









A 

M, 



,( t) 1 



( 4. 27) 



A formula for the conditional variance of M^(t) 
original sample data is derived in Appendix 



expression is given by a 

Var[ M.( t ) ] = i-x? - [ ” ^ 






-2t[ ----] [ --SLx-l } 



A A 
+• 2 r 1 



given the 
C. This 



( 4. 28) 



Notice that the expected value of the conditional variance 
of the bootstrap estimate is approximately equal to the 
variance of the nonparametric estimate of M^j(t) which is 
derived in Appendix A. 

So far we have considered the problem, where the arrival 
rate is known. The bootstarp methodology also applies if 
the arrival rate is unknown and is estimated from 
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interarrival data. Suppose the data consist of a random 
sample X=(X, ,X_j_ , . . . ,X^ ) from unknown service time 

distribution F and an independent sample Y=( Y, ,Y^ ,Y^) 

from the exponential interarrival time distrbution G with 
unknown parameter . One bootstrap procedure to estimate 
the expected number of customers being served at time t is 
to construct and , the empirical probability 

distribution corresponding to F and G. Bootstrap samples 
X-^~ F^ , i=l,2,...n, Y;'*~ Gp^, j=l,2,...n, are independently 

drawn, an estimate of Mg(t) 

mJ( t) = ? A + -Mn - i I(xf<t)] t] (4.29) 

is calculated. As before there are a large number B of 
bootstrap replications. For this case, the bootstrap 
estimate of M^(t) and its variance are still given by 
equations 4. 26 and 4. 27. There appears to be no closed form 

A. 

of the analytical variance of (t) in this case. Now we 
will describe methods to obtain approximate confidence 

A 

intervals for the bootstrap estimate M^(t). 

1. The Percenti le Method 

A simple method for assigning approximate confidence 
intervals to the nonparametric estimate M^j(t) is as follows: 
Let 

C(t) = (4.30) 

B 

be the cumulative distribution function of the bootstrap 
distribution of M,.j(t); B is the number of bootstrap 
replications. For a given 0<c<<0.5, define 
UoO = c'(oO , U(c^) = c^d-oo 

A. 

Usually denoted simply by L and U. This definition runs 
into complications when we actually try to compute quantiles 

A A 

L and U from a set of bootstrap replications. To overcome 
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these difficulties, we order the bootstrap replications from 

jH,*) 

smallest to largest, obtaining the sorted data (t), for 
i=l to B. Letting represent any fraction between 0 and 
1; take Q(c<) to be (t) whenever Q is one of the 

functions , for i=l to B. Thus L(o<) turns out to 

be the (B and U(o(,) to be the (B * ( l-“<) 

+ 0.5)^'^ M’^*'^(t). The percentile method consists of taking 

[ L(o<) , U(o0 ] (4.31) 



as an approximate l-2o( confidence interval for t) since 

A ^ 

<X=C(L), l-o<=C(U), the percentile method interval consists 

of the central 1-2 o< proportion of the bootstrap 
distribution. 

2. The Bias- corrected Percentile Method 

Efron(1980) suggests the following bias correction 
for the percentile confidence interval procedure [Ref. 11]. 

A 

He argues that if Mq( t) is not the median of the bootstrap 
replication distribution, then a bias correction to the 
percentile method is called for. To be specific, define 



z o = i‘'[C(M^(t))j (4.32) 

where C(t)= ^ as in equation 4.30, and is the 

cumulative distribution function for a standard normal 
variate. The bias corrected percentile method consists of 
taking 

[C''[£(2zo- z^)} , c'[£(2z^+ z^)}] (4.33) 

as an approximate 1-2 o( central confidence interval for 

A 

Mo^(t). Here Zo<^is the upper point for a standard normal 
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Notice that if Mg(t) is the median of the bootstrap 
distribution then Zq=0 and equation 4. 33 reduces to equation 
4.31, the uncorrected percentile interval. However, even 
small differences of Pr(Mj^(t)f, M^( t) j from 0.5 can make 
equation 4. 33 much different from equation 4. 31. 
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V. SIMULATION RESULTS 



The purpose of the simulation in this chapter is to 
assess the performance of the nonparametic estimation 
methods, the jackknife and the bootstrap. Since the 
estimate of M(t), the mean number of customers being served 
at time t, is a function of the customer arrival rate and 
the integral of the survivor function of the service time 
distribution, two simulations cases are done. The first 
simulation case was performed to estimate M^j(t), the 
nonparametric estimate of M( t) , as a function of the service 
times with the arrival rate assumed to be known and set 
equal to 1. For this case, the jackknife and bootstrap 
estimate of the variance were derived in the chapter IV, and 
compared with the numerical estimate obtained by the 
simulation. The second case assumed that the customer 
arrival rate is also unknown and must be estimated using 
interarrival times. 

In each replication of the simulation for case 1, 50 
independent service times from a specified service time 
distribution were generated. For the bootstrap procedure, 
500 bootstrap replications were performed. The simulation 
was replicated 300 times. For the purposes of comparison, 
we considered four types of service time distributions, 
which were the exponential, the mixed exponential, the 
gamma, and the lognormal distribution. The arrival process 
is known to be Poisson process with known rate X =1. The 
same generated service times were used for each estimation 
procedure in a replication. This reduces the variability of 
the differences in performance between the procedures. All 
programming was done on IBM 3033 computer at the Naval 
Postgraduate school using the LLRANDOMII, random number 
generating package [ Ref. 6] . 
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TABLE III 



STATISTICAL DATA OF ESTIMATE OF M(T) 
N=50 R=300 (B=500) 





True 

value 


Parametric 


Nonparametric 


Correct 


Errors 


Jack 


Boot 


Exponential 
( M = 2) 


0. 7869 


0. 7830 
( 0. 0268) 


- 


0. 7820 
( 0. 0451) 


0. 7819 
( 0. 0452) 


Mixed expon. 

— 2,Wi:— . V5,p, — . ^ 


0. 5992 

2) 


0. 5871 
( 0. 0023) 


0. 6246 
( 0. 0026) 


0. 5991 
(0. 0512) 


0. 5989 
( 0. 0512) 


II O 
II 

tSJ 


0. 8963 


0. 8952 
( 0. 0009) 


0. 7861 
( 0. 0010) 


0. 8965 
( 0. 0325) 


0. 8967 
( 0. 0325) 


Lognormal 
(§=. 193, <r-=l) 


0. 8094 


0. 8140 
( 0. 0021) 


0. 7844 
( 0. 0020) 


0. 8139 
( 0. 0383) 


0. 8138 
( 0. 0383) 



Table III presents the results of several estimation 
methods when the arrival rate is given and equal to 1. The 
top of each cell gives the mean estimate of M( t) at time 
t=l, where M( t)= 0F( s)ds. The bottom part of each cell 
gives the standard deviation of the estimate. For service 
time distributions other than exponential, a parametric 
estimate based on an erroneous exponential model is also 
given. The estimate in the case of an exponential model is 
[ 1-EXP { -t/u 1 ] where aa is the mean service time. For each 
service time distribution, the standard deviation of the 
parametric estimate of M(t) is smaller than that of the 
nonparametric estimate of M(t). That is, the efficiency of 
the parametric estimation method is better than the 
efficiency of the nonparametric estimation method. However, 
the results of a parametric fit assuming an erroneous 
exponential model show the worst performance. The true 
value of M(t) is not included within plus or minus three 
standard deviations of the erroneous estimate M(t). In the 
table, the nonparametric estimation methods seem to perform 
well in all cases with the cost of an inflation of variance. 
Hence the nonparametric estimation method is to be preferred 
when the service distribution is unknown. 
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To illustrate the efficiency of the nonparametric 
estimation methods, we simulated two possible ways to 
construct the approximate confidence interval for M^i(t) for 
the bootstrap and the jackknife methods. Those ways are 
presented in chapter IV. For the jackknife estimation 
method, one procedure was to construct the confidence 
interval with the regular degrees of freedom, n-1, and the 
other used the reduced degrees of freedom, which is the 
number of different pseudo values. For the bootstrap 
estimation method, one way used the percentile method by the 
Monte Carlo process, and the other used the bias-corrected 
percentile method; there were 500 bootstrap replications. 
Nominal 68%, 80%, and 90% confidence intervals were 
constructed for each replication using each method. It was 
noted whether the confidence interval formed by a given 
method covered the true value M(t). The entire process was 
independently replicated with R=300 times. From these R 
replications we computed, for each method, the proportion p 
of the R confidence intervals which contained M(t), as well 
as the average length of the confidence intervals. If a 

A 

method was performing adequately, p should be near l-c>< , and 
a small mean length is desirable. 

Tables IV to VII show the simulation results of several 
confidence interval procedures for four types of service 
time distribution; the exponential, the mixed exponential, 
the gamma, and the lognormal. The arrival process is 
Poisson with known arrival rate \=1. 

In order to compare the performance of these procedures 
to the normal confidence interval procedure, simulations 
were conducted, and nominal 68%, 80%, and 90% confidence 
limits were constructed for time t=l for each replication. 
The normal confidence interval procedure is based on the 
order statistics of the service times. By the central limit 
theorem, the distribution of M^(t) is asymptotically normal 
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TABLE IV 



COVERAGE AND LENGTH OF 100( 1-'^)% C. I 
FOR EXPONENTIAL WITHaa=2, \ =1 ATT=1 







68 ^ 


/ 


80 % 




90 % 








Length 

(s.d) 


C. R. 


Length 

(s.d) 


C. R.. 


Length 

(s.d) 


C. R. 


Normal C. I. 


0. 0888 


19. 33 
69. 33 


0. 1147 


12. 00 
83. 67 


0. 1477 


4. 67 
92. 00 


Procedure 


( 0. 0089) 


11. 33 


( 0. 0101) 


4. 33 


( 0. 0141) 


3. 30 




Reduced 


0. 0911 


17. 67 
71. 33 


0. 1187 


11. 00 
85. 33 


0. 1548 


3. 33 
94. 00 


Jack- 


d. f 


( 0. 0089) 


11. 00 


( 0. 0102) 


3. 67 


( 0. 0142) 


2. 67 


knife 


Regular 


0. 0897 


18. 38 
70. 67 


0. 1163 


11. 33 
84. 33 


0. 1505 


4. 00 
93. 00 




d. f 


( 0. 0090) 


11. 00 


( 0. 0103 ) 


4. 33 


( 0. 0144) 


3. 00 




Percen- 

tile 


0. 0884 


18. 67 
69. 00 


0. 1132 


12. 00 
82. 00 


0. 1462 


4. 33 
92. 00 


Boot- 


method 


( 0. 0098) 


12. 33 


( 0. 0112) 


6. 00 


( 0. 0154) 


4. 67 


strap 


Bias- 

correct 


0. 0887 


16. 67 
71. 00 


0. 1137 


10. 67 
83. 00 


0. 1467 


2. 67 
92. 67 




method 


( 0. 0098) 


12. 33 


( 0. 0112) 


6. 33 


( 0. 0154) 


4. 67 



distributed as the number of data points n -^oo . Thus, the 
100( l-c<)% normal confidence interval is given by 



M^(t) 


+ 


z,.ayVar[M^( t)] 


(5. 1) 


where z 


'-1 


is the upper 1-^ point of the 


standard normal 



distribution and Var[M>j(t)] is given by equation A. 14 in 
appendix A. 

Each cell in the tables contain the average and standard 
deviation of confidence interval length; and the proportion 
of intervals that are too high, ( e. g. M (t)<L), where L is 
the lower bound of interval; the proportion of intervals 
covering the true value M(t), p; the proportion of interval 
that are too low, ( e. g. M(t)>U), where U is the upper bound 
of interval. Table IV is for the exponential service time 
case with with H=2; Table V is for the mixed exponential 
service time case with Mi=2, 75, and P, =0.2; Table VI 
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TABLE V 



COVERAGE AND LENGTH OF 100( l-o<)% C. I. FOR 
MIXED EXPONENTIAL WITH w, =2, Mx=. 75, P, =. 2, \ =1 AT T=1 







68 ? 




80 5 




90 ^ 








Length 

(s.a) 


C. R. 


Length 
( s. d) 


C. R. 


Length 
( s. d) 


C. R. 


Normal C. I. 


0. 0914 


14. 67 
71. 67 


0. 1169 


10. 00 
83. 00 


0. 1509 


5. 00 
90. 67 


Procedure 


( 0. 0084) 


13. 67 


( 0. 0102) 


7. 00 


( 0. 0143) 


4. 33 




Reduced 


0. 0936 


14. 00 
73. 00 


0. 1208 


9. 67 
83. 33 


0. 1581 


4. 33 
91. 67 


Jack- 


d. f 


( 0. 0085) 


13. 00 


( 0. 0103) 


7. 00 


( 0. 0143) 


4. 00 


knife 


Regular 


0. 0923 


14. 67 
72. 00 


0. 1185 


10. 00 
83. 00 


0. 1538 


4. 33 
91. 33 




d. f 


( 0. 0085) 


13. 33 


( 0. 0103) 


7. 00 


( 0. 0145) 


4. 33 




Percen- 

tile 


0. 0911 


14. 33 
71. 67 


0. 1156 


11. 33 
81. 33 


0. 1490 


4. 33 
89. 67 


Boot- 


method 


( 0. 0094) 


14. 00 


(0. 0112) 


7. 33 


( 0. 0153) 


4. 33 


strap 


Bias- 

correct 


0. 0913 


14. 00 
71. 00 


0. 1161 


9. 00 
83. 00 


0. 1495 


4. 00 
89. 67 




method 


( 0. 0094) 


15. 00 


( 0. 0112) 


8. 00 


( 0. 0153 ) 


6. 33 



is for the gamma service time case with 3=1 and K=2; and 
Table VII is for the lognormal service time case with 
§ =0. 193 and <T'^=1. 

The overall examination of the tabulations of confidence 
limit coverage and also the average and standard deviation 
of confidence interval length suggest that the bootstrap 
procedure is slightly better than the jackknife procedure; 
however, the difference is negligible. The normal 
confidence interval is also about the same as the jackknife 
and bootstrap procedures indicating that a sample size of 50 
is large enough for the central limit theorem approximation 
to be adequate. All procedures produce almost the same 
average length of confidence interval with a good coverage 
rate, which falls within ± 2 l-o< . Although the 

1 2oO 

method of the reduced degrees of freedom used in the 
jackknife and the bias-correct percentile method applied in 
the bootstrap improved the coverage rate, the variance was 
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TABLE VI 



COVERAGE AND LENGTH OF 100(l-o^)% C. I. 
FOR GAMMA WITH ^=1, K=2 , \=1 AT T=1 







68 ? 




80 ? 


/ 


90 5! 








Length 
( s.d) 


C. R. 


Length 

(s.d) 


C. R. 


Length 
( s.d) 


C. R. 


Normal C. I. 


0. 0595 


22. 00 
68. 33 


0. 0774 


12. 33 
80. 33 


0. 0989 


10. 67 
86. 33 


Procedure 


( 0. 0100) 


9. 67 


( 0. 0120) 


7. 33 


( 0. 0174) 


3. 00 




Reduced 


0. 0617 


21. 00 

70. 00 


0. 0813 


11. 33 
82. 67 


0. 1062 


8. 33 
89. 00 


Jack- 


d. f 


( 0. 0102) 


9. 00 


( 0. 0122) 


6. 00 


( 0. 0178) 


2. 67 


knife 


Regular 


0. 0601 


21. 67 
69. 33 


0. 0785 


12. 33 
80. 67 


0. 1008 


10. 33 
86. 67 




d. f 


( 0. 0102) 


9. 00 


( 0. 0121) 


7. 00 


(0. 0177) 


4. 00 




Percen- 

tile 


0. 0589 


21. 33 
69. 97 


0. 0766 


12. 00 
80. 33 


0. 0981 


9. 33 
86. 67 


Boot- 


method 


( 0. 0091) 


9. 00 


( 0. 0122) 


7. 67 


( 0. 0179) 


3.00 


strap 


Bias- 

correct 


0. 0595 


19. 33 
69. 33 


0. 0777 


10. 67 
81. 00 


0. 0990 


8. 67 
87. 00 




method 


( 0. 0100) 


11. 33 


( 0. 0121 ) 


8. 33 


( 0. 0180) 


4. 33 



inflated. Furthermore, the amount of improvement was small 
and not significant. Hence, the original procedures for 
constructing the confidence interval for the jackknife and 
bootstrap are preferred in this case. Note that the 
coverage rates are skewed left slightly but almost balanced. 
It is a reason that the normal confidence interval procedure 
performs well. 

Results will now be reported for the simulation of the 
case in which the arrival rate of the Poisson process is 
also unknown and must be estimated from interarrival time 
data. More computations are required for this case; 
however, the procedure is same. Each replication of the 
simulation generated 50 independent service times and 50 
independent exponential interarrival times having mean 1. 
Confidence intervals were computed using both separated and 
paired jackknife procedures and the percentile method for 
the bootstrap. The number of bootstrap replications was 
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TABLE VII 



COVERAGE AND LENGTH OF 100(1-®^)% C. I. 

FOR LOGNORMAL WITH §=.193, \ =1 AT T=1 







68 % 




80 % 




90 % 








Length 
( s. d) 


C. R. 


Length 
( s.d) 


C. R. 


Length 
( s. d) 


C. R. 


Normal C. I. 


0. 0769 


16. 67 
71. 00 


0. 1007 


10. 00 
78. 00 


0. 1272 


6. 67 
89. 33 


Procedure 


( 0. 0075) 


12. 33 


( 0. 0096) 


12. 00 


( 0. 0126) 


4. 00 




Reduced 


0. 0787 


16. 33 
71. 33 


0. 1208 


9. 67 
79. 67 


0. 1330 


6. 33 
89. 67 


Jack- 


d. f 


( 0. 0076) 


12. 33 


( 0. 0097) 


10. 67 


( 0. 0127) 


4. 00 


knife 


Regular 


0. 0763 


16. 67 
70. 33 


0. 1021 


10. 00 
79. 00 


0. 1297 


6. 67 
89. 33 




d. f 


( 0. 0076) 


13. 00 


( 0. 0097) 


11. 00 


( 0. 0128) 


4. 00 




Percen- 

tile 


0. 0763 


16. 67 
69. 33 


0. 0998 


10. 33 
77. 00 


0. 1258 


7. 33 
88. 67 


Boot- 


method 


( 0. 0084) 


14. 00 


( 0. 0105) 


12. 67 


( 0. 0133) 


5. 00 


strap 


Bias- 

correct 


0. 0768 


16. 00 
68. 33 


0. 1004 


8. 67 
78. 33 


0. 1263 


6. 33 
88. 67 




method 


( 0. 0084) 


15. 67 


( 0. 0106) 


13. 00 


( 0. 0133) 


5. 00 



1000. Nominal 68%, 80%, and 90% confidence limits were 
computed for each replication. The simulation was 
replicated 300 times. 

Tables VIII to X report the results of the simulation. 
The quantities in the left part of each cell are the average 
and standard deviation (within parenthesis) of coverage 
interval length. The right part of each cell contains three 
quantities; the top value is the proportion of intervals 
that are too high; the center value is the proportion of 
intervals that cover the true value, p; and the bottom part 
is the proportion of intervals that are too low. 

In Table VIII (the case of 68% C. I. ), the average length 
from the bootstrap shows outstanding performance with a 
small value of standard deviation. The paired jackknife 
procedure performs as well as the bootstrap procedure. This 
procedure reduced the standard deviation by more than half 
of that in the separated jackknife procedure, and also 
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improved the coverage rate. From the results of coverage 
rate in the table, it can be recognized immediately that the 
jackknife estimate, regardless of the application method, is 
often too low, while the bootstrap estimate tends may a 
little too high but is almost balanced in the number of 
confidence intervals that are too high or too low. It is 
the reason that the bias-corrected percentile method was not 
required in this case. 

In this simulation, all the coverage rates fall within 
± 2 of 1” / ( 62.61, 73.38). Note that the average 
length of the confidence interval in the gamma service time 
case is the highest. When the arrival rate was known, the 
gamma service time case had the smallest average length. 
This indicates that the variability of the estimated arrival 
rate may be the dominate effect in the width of the 
confidence interval. 

The results of Tables IX and X support the facts of 
discussion about Table VIII, This is the reasonable since 
the same random numbers were used to compute these 
confidence interval. The presentation of Table IX is 
exactly the same as the case of table VIII. All the 
coverage rate are fall within (75.38, 84.61), though the 
value of coverage rates fluctuate over the service time 
distribution cases. Obviously the paired jackknife 
procedure performs very well. The bootstrap procedure still 
has the best performance; however the value of coverage rate 
fluctuates greatly for the different service time 
distributions. For the 90% confidence interval case 
reported in Table X, some coverage rate fall outside the 
range (86.53, 93.46). The separated jackknife procedure 
produces low coverage rates outside of (86.53, 93.46), 
except for the exponential service time distribution. The 
paired jackknife procedure improved the coverage rate 
tremendously. Although the average lengths in the paired 
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jackknife procedure are slightly bigger than those in the 
bootstrap procedure, the overall performance is better than 
the bootstrap. Furthermore, the procedure in the case of 
gamma service time distribution, the bootstrap produced one 
coverage rate outside of (86.53, 93.46). However, this 
could be due to sampling fluctuation. 

In general, all the confidence interval procedures 
performed very well for the exponential service time case, 
regardless of the level of the confidence interval. The 
procedures also worked well in general to produce 68% and 
80% confidence interval. However, performance was more 
variable in the 90% confidence interval case. In most 
cases, the average length produced by the bootstrap 
procedure is the smallest, but the value of the coverage 
rate fluctuates for different service time distribution. 
The overall examination of the tabulations suggests that the 
paired jackknife procedure performs very well compared to 
the separated jackknife procedure and in some cases shows 
better performance than the bootstrap procedure. 
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TABLE VIII 



COVERAGE AND LENGTH OF 68% C. I. 

WITH UNKNOWN ARRIVAL RATE (N=50, R=300, B=1000) 







Jackknife 




Bootstrap 




Separated 


Paired 


Percent! le 




Length 
( s. d) 


C. R. 


Length 

(s.d) 


C. R. 


Length 
( s.d) 


C. R. 


Exponential 


0. 2674 


7. 00 
70. 00 


0. 2525 


9. 67 
70. 33 


0. 2436 


18. 33 
66. 67 


( AA= 2) 


( 0. 1282) 


23. 00 


( 0. 0570) 


20. 00 


( 0. 0510) 


15. 00 


Mixed exponen 
( Ml =2, Mi-=. 75 


0. 1999 


12. 00 
65. 67 


0. 2077 


15. 00 
68. 00 


0. 2009 


21. 67 
64. 00 


( 0. 0832) 


22. 33 


( 0. 0429) 


17. 00 


( 0. 0398) 


14. 33 


Gamma 


0. 2917 


10. 00 
68. 00 


0. 2746 


13. 33 
67. 67 


0. 2632 


21. 67 
66. 00 


( /3=1, k=2) 


( 0. 1376) 


22. 00 


( 0. 0599) 


19. 00 


( 0. 0557) 


12. 33 


Lognormal 


0. 2533 


7. 67 
72. 00 


0. 2513 


10. 00 
: 72. 33 


0. 2415 


19. 33 
69. 67 


(§=. 193, cf^=l) 


( 0. 0998) 


20. 33 


( 0. 0486) 


17. 67 


( 0. 0461) 


11. 00 



TABLE IX 

COVERAGE AND LENGTH OF 80% C. I. 

WITH UNKNOWN ARRIVAL RATE (N=50, R=300, B=1000) 







Jackknife 




Bootstrap 




Separated 


Paired 


Percentile 




Length 
( s. d) 


C. R. 


Length 
( s. d) 


C. R. 


Length 
( s.d) 


C. R. 


Exponential 


0. 3447 


5. 33 
79. 67 


0. 3282 


8. 00 

82. 00 


0. 3162 


12. 67 
83. 00 


( M = 2) 


( 0. 1329) 


15. 00 


( 0. 0641) 


10. 00 


( 0. 0585) 


4. 33 


Mixed exponen 
( M. =2 , M^=. 75 


0. 2558 


75. 67 


0. 2688 


5. 67 
82. 33 


0. 2553 


12. 67 
83. 00 


P, =-^) 


( 0. 1088) 




( 0. 0569) 


12. 00 


( 0. 0490) 


4. 33 


Gamma 


0. 3646 


5. 00 
78. 67 


0. 3504 


6. 67 
81. 67 


0. 3375 


16. 33 
75. 67 


( ^=1, k=2) 


( 0. 1509) 


19. 33 


( 0. 0703 ) 


11. 67 


( 0. 0657) 


8. 00 


Lognormal 


0. 3259 


4. 67 
77. 67 


0. 3231 


6. 33 
79. 33 


0. 3115 


14. 67 
75. 67 


(§=. 193, <T =1) 


( 0. 1279) 


16. 67 


( 0. 0624) 


14. 33 


( 0. 0585) 


10. 00 
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TABLE X 



COVERAGE AND LENGTH OF 90% C. I. 

WITH UNKNOWN ARRIVAL RATE (N=50, R=300, B=1000) 







Jackknife 




Bootstrap 




Separated 


Paired 


Percentile 




Length 
( s. d) 


C. R. 


Length 
( s. d) 


C. R. 


Length 
(s. d) 


C. R. 


Exponential 


0. 4464 


1. 33 
90. 00 


0. 4231 


2. 67 
92. 33 


0. 4091 


8. 67 
89. 00 


( M = 2) 


( 0. 1680) 


8. 67 


( 0. 0844) 


5. 00 


( 0. 0761) 


2. 33 


Mixed exponen 
(M, =2. 75 


0. 3199 


0. 67 
83. 33 


0. 3387 


2. 33 
88. 00 


0. 3301 


8. 00 
86. 67 


P, =. ^) 


( 0. 1213) 


16. 00 


( 0. 0679) 


9. 67 


( 0. 0609) 


5. 33 


Gamma 


0. 4891 


2. 00 
85. 33 


0. 4586 


3. 33 
88. 67 


0. 4416 


9. 00 
85. 67 


( /S=l, k=2) 


( 0. 2321) 


12. 67 


( 0. 1053) 


8. 00 


( 0. 0959) 


5. 33 


Lognormal 


0. 4371 


1. 67 
85. 33 


0. 4228 


3. 00 
90. 00 


0. 4074 


7. 00 
88. 67 


(§=. 193, (f=l) 


( 0. 1974) 


13. 33 


( 0. 0920) 


7. 00 


( 0. 0847) 


4. 33 
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VI. SUMMARY AND CONCLUSIONS 



This thesis consider the problem of estimating M(t), the 
mean number of customers being served at time t for an 
M/G/oo queue, using service time and interarrival time data. 
It is assumed that there are no customers being served at 
time 0. Two cases are considered. In one the parametric 
form of the service time distribution is assumed known. In 
this case M( t) is a function of the estimated parameters. 
In the situation in which the arrival rate of the Poisson 
process is also assumed known and the parametric form of the 
service time distribution is exponential, approximation to 
the bias and variance of the estimate are derived. Further, 
simulation is used to study a normal confidence interval 
procedure. 

For the other case the parametric form of the service 
time distribution is unknown. The empirical distribution of 
the service time distribution is used in the estimate of 
M(t). In the situation in which the arrival rate \ is 
assumed known, the distribution of the estimate is derived 
in Appendix A. The bootstrap and jackknife estimates with 
known are studied in Appendix B and C. Simulation was used 
to assess the performance of confidence interval procedures 
using a normal approximation the jackknife and the 
bootstrap. The simulation results for the case in which the 
arrival rate \ is known indicate that: 

(1) The parametric estimation method appears the most 
powful method when the parametric assumption is 
correct, but the performance is seriously degraded if 
the assumption is not appropriate. 

(2) When an erroneous parametric (exponential) model is 
assumed, the initial estimates of mean number in 
service are poor. However, as t -^oo , the erroneous 
parametric estimate approaches the same value as the 
other estimates. This is because as t -»oo all the 
estimates approach the sample mean of the service 
time data. 
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( 3 ) 



The estimate obtained by using the empirical 
distribution is unbiased with a larger variance than 
a parametric estimate based on a correct model. 

(4) Simulation results indicate there is not much 
difference between jackknife and bootstrap confidence 
interval procedures. 

(5) The nonparametric normal confidence interval 

procedure performs as well as the procedure in ( 4) 
since the distribution of the estimate is almost 
symmetric. The improvement by the use of adjusted 
degrees of freedom in the jackknife ana the 

bias-corrected percentile in the bootstrap is small. 

We now discuss the simulation results for the case in 
which the arrival rate for the Poisson process is also 

unknown and is estimated using interarrival time data. The 
service times are generated from four types of distribution. 
The percentile method for the bootstrap and paired and 
separated techniques for the jackknife were used to 
construct the confidence intervals. Tables I and II, which 
are the results of a parametric confidence interval 

procedure in chapter III, are compared with the results of 
the nonparametric confidence interval procedures. The 
simulation results indicate that: 

(1) The nonparametric confidence interval procedure works 

as well as the parametric case, even though the 

length of the confidence interval is wider than the 
parametric one. 

(2) In the overall examination, the percentile method of 
the bootstrap shows the best performance. The paired 
jackknife procedure also has similiar results to the 
bootstrap approach. The results of these two 
nonparametric procedures show the almost same level 
of performance with the parametric one. However, the 
separated jackknife procedure produces poor results. 

(3) The results of the jackknife procedures produce 
intervals that are always biased upward. Efron 
(1981) reported similiar results. [Ref. 13] 

(4) Since the bootstrap procedures require a large amount 

of computation, the jackknife is the method of choice 
if one does not want to do the bootstrap 

computations. 

In general, the nonparametric methods of the bootstrap 
and the jackknife performs very well, regardless the 
complexing of the estimation problem. Of cource, if the 
parametric estimation method can be applied, the results are 
clearly superior. However, the application of the 
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parametric estimation is a highly limited because the 
parametric assumption is often difficult to verify. When 
the estimate is simple enough, which is the nonparametric 
estimate when the arrival rate is known, and the asymptotic 
distribution of estimate can be obtained, the nonparametric 
normal confidence interval procedure performs well, and more 
complicate computations such as the jackknife and the 
bootstrap method are not required. However, the jackknife 
and the bootstrap method have a good performance for the 
more complicated problem in which the arrival rate is 
unknown. The bootstrap confidence intervals show the best 
performance but the paired jackknife procedure achieve the 
same level of performance with less computation than the 
bootstrap in this problem. 
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APPENDIX A 

CALCULATING THE BIAS AND THE VARIANCE OF M^(T) WITH KNOWN 

ARRIVAL RATE 



It will be assume that the rate of Poisson process \ is 
known and is equal to 1. The nonparametric estimate Mvi(t) 
is given by 

Mvi(t) = 5^^ [ 1-F( s)]ds (A. 1) 



Using the empirical cumulative 
nonparametric estimate is 



t) 



--Xs, .+ 

^ (a) 



t 



density function Fr^ , the 



(A. 2) 



where the observation s are the order statistics of the 
indeprndent and identically distributed service times with 
unknown distribution F. To find the distribution of Mt 4 (t), 
we will study the distribution of 

Let be independent, identically distributed 
random variable with distribution function F, Let N be the 
number of X;'s which are less than t. Let X^^^ denote the i"^*^ 
smallest X;. By the definition of conditional probability. 



x|N^=l} 






(A. 3) 



for x<t. Since the random variable N-^ 
distribution with a parameter F(t), we 
equation A. 3 to obtain 

Pit') C 



nfL 

Fit) 



has a binomial 
can rewrite the 



(A.4) 
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Let f(x) be the density of the distribution function F. The 
conditional probability given N^=2, 

P { X , ^ , + dx, , x^^ x^+ dx^ I Nt =2 ) 



_ 2 

"y'o tW 



(A. 5) 



for X ,<x_j^. 

A 

Given N^=K, the conditional distribution of the values 
of the unordered X; that lie in (0,t] is that of independent 
random variables with distribution function , for 0^x<t. 

Thus, given N^=K, M^t) has the same distribution as a 
constant plus the sum of K independent identically 
distributed random variables. Thus the expectation of M^(T) 
can computed by the property of conditional expectation. 



E[M^i(t)] = E[E[M^j( t) |Nt] ] 
Given N-t=K/ 



(A. 6) 






(A.7) 



Since the random variable has a binomial distribution 

with the parameter F(t), 



= s. 



xF( dx) + [ 1-F( t ) ] t 



(A. 8) 



To check the bias of the estimate (t), using the 

integration by part of equation A. 1, the true estimate is 
given by 



M^(t) = t[l-F(t)] + £^tF(dx) 



(A.9) 
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Using conditional 



Thus the estimate M^j(t) is unbiased, 
expectations, the variance is computed by 
Var[M^i(t)] = E[(M^(t) - E[Mg(t)])2] 



= E[(M^(t) - E[M^(t) |Ntl )2] 



+ 2E[(M^(t) - E[M^t) |N.tn(E[M^(t) IN^] - E[M^(t)])] 
+ E[ (E[M^(t) |N.tl -E[M„(t)])2] (A.IO) 



Computing each of the individual terms, we obtain 
E[ (My( t) - E[M^( t) 1%] )2] 



F_^t) 




. r 



X 



rosx) 



(A. 11) 



and 

E[(E[M^i(t) IN^l - E[MJt)] )2] 

= -- F(t)F(t)[ Tx - t] 2 (A. 12) 

where F is the survival distribution function of F. The 
second term of right term of equation A. 10 turn out to be 
zero. Thus the variance of Mfi(t) is obtained by sum of two 
equations. The resulting variance is 

Var[M^j(t)] = { [^^x2F(dx) - [ C xF( dx) ] 2 + t2f( t )F( t ) 

- 2tF( t)[ t^xF(dx)] ] (A. 13) 



Notice that the variance estimate of Mjj(t) is equal to 
"Var(X) as t->oo . A nonparametric estimate of Var[M»j(t)] is 



Var[M^(t)] 






(.i) 



- 



V. 



r\ 



[ 1 - — ■ 
n • 
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where K is the number of service times that are less than t. 



APPENDIX B 

JACKKNIFE ESTIMATE OF NONPARAMETRIC ESTIMATE 



It will be assumed that the arrival rate \ is known and 
equal to 1. The nonparametric estimate is given by 



M. 



A 

K 



A 

n- K 



(t) = --Xs,, + 

' ^ y\ C*> r\ 



(B. 1) 



where s are the order statistics of the service times and 

A 

assumed that the random variable K exists such that 

A A 

The estimate M ( t) for the data set 



<J^) ((«■» 



obtained by delating the i'^'' point from the sample is 
' ^ -) - 1< 






— s ^ + 

-J=\ Cj> 



r\-\ 

A 

n- K 



--- 4. s^5 + t 

i\-\ n-\ 

0’ ' 



if i > K 



if i £ K 



(B.2) 



The pseudo-value M;(t) is computed by 



M-(t) = nM^^j(t) - (n-l)M^^(t) 



(B.3) 



where M^^^(t) is the estimate of M,j( t) based on all the data. 



substituting the estimate M^(t) in equation B.3, we get 



M-(t) = 



1 ^Ca') 



if i ^ K 



t if i > K 

Since the jackknife estimate is the average 
pseudo-values, the estimate is given by 



^ A 

M ^ ^ n-K , 

M-ft) = --Xs,.. + t 



(B.4) 



of the 



(B. 5) 
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Thus, the jackknife estimate is the same as the original 
nonparametric estimate of Mfj(t). The jackknife estimate of 
variance is 



Var[M3(t)] = .L.[ J- Sl{M^(t)} 



2 _ 






(B. 6) 



where 
and A 

Thus the variance cai\ be rewri ^cen as -\ ^ 

Var[M3(t)]= 






A 

n" ^ 



(B.7) 



Comparing equation B. 7 with equation A. 14, it is seen that 
the jackknife estimate of variance is greater than the 
estimate of equation A. 14 by a multiplicative constant 
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APPENDIX C 

CONDITIONAL DISTRIBUTION OF BOOTSTRAP ESTIMATE 



It will be assumed that the arrival rate \ is known, 
and is equal to 1. Let S,,S^,..,Sy^ denote random service 
times from a CDF F, and let S,.^S <...^S ^t^ S ^..£S,_^ denote 
the corresponding order statistics. The nonparametric 
estimate of j^F(s)ds is 



M«( t) 



I ^ n-K 

— — + — — — 

^ IT| 



t 



(C. 1) 



Let B^, i=l to n, be independent random variables having the 
same distribution as draws with replacement from ( s, , Sj^ , . . , s^ 
), and let / i = l to n, be the corresponding order 
statistics. A bootstrap realization of the nonparametric 
estimate is 

M^(t) = ^ -;-ln-|l(l;,,it)lt (C.2) 

where 

I(x^t)=^l ifxlt 

I 0 otherwise 

A 

To compute the distribution of M^(t), the Laplace transform 

A 

is used [Ref. 12]. The Laplace transform of Mg(t) is 

E[ exp( -f M^(T) )] = E[E[exp(-5Ma(t))lN^,j] (C.3) 



by the property of conditional expectations, where N^t is the 
number of bootstrap samples which are less than or equal to 
t. We compute the right hand side of equation C. 3 
separately. First 

E[exp(-§ Mg(t)) lN<-= 1] 
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(C.4) 



= exp( t ) [ exp( - ^ exp( - § ) ] 



is computed, where the random variable Nt has a binomial 

A 

distribution with the parameter F(t)= . Thus, from 

A 

equation C. 3 the Laplace tra^nsform of M^( t) is 



E[exp(- M (t)] = exp( t) exp( -1-^-)^ -I^exp( 






= exp(-|t)[-'^exp(|-^)^-^exp(-j-^^^) + (C. 5) 



Let us define a random variable Y having the following 
distribution 



Y = 



S(/:> 1 . ^ 

w. p — If 1 < K 



(C.6) 



if i > K 



^ w. p __ 

Then the Laplace transform of Y is 

k 



E[exp(-fY)] = ” Xexp( ) + -'^-exp( ) (C. 7) 



Thus, (t) has the same distribution as the sum of n 
independent random variables having the same distribution as 
Y. For the fixed time t, given the order statistics, 

<. . <S <t<S <. . <S , the expectation of A^(t) is written by 
E[ Mg( t ) I data] = nE[Y|data] 



A 

I J< n-t 

= -- ^s +-jr-t 

n 



(C.8) 



Thus, the bootstrap estimate is asymptotically an unbiased 
estimate of M^j(t). The variance is 



Var[ M„( t ) I data] = nVar[ Y] 



(C.9) 



The variance of Y can be derived using the equation C. 6 
Since 
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E[ Y2] 



( C. 10) 



J_ X [ ---i ] 2 + JL ) 2 
n ^ n n ^ ^ ^ 



Hence the asymptotic bootstrap variance estimate of MN(t) is 
given by 



Var[ Mq( t ) I data] = -'-{--^s? - [-J-Xa.]^ + t^-^-- --- 
B' ' ' n n x-i h n r\ 






i ^ 



K n-K 



r\-i< 



A 



- 2t [-;,-zs ]2} 



(C. 11) 



That is, 
the same 



the asymptotic bootstrap estimate of variance 
as the nonparametric estimate (equation A. 14). 



is 



70 



LIST OF REFERENCES 



1. Gross, D. and Harris, C. M. , Fundamenta 1 s of Queueing 
Theory . John Wiley and Sons, 1974. 



2. Brown, M. and Ross, S. M. , "Some Results for Infinite 
Server Queues", J. Appl . Prob . . 6, pp. 604-611, 1969. 



3. Takacs, L. , Introduction to the Theory of Queues . 
Oxford University Press, T962T 



4. Cox, D. R. < "Some Problems of Statistical Analysis 
Connected with Congestion , Proceedings of the 
Symposium on Congestion Theory , pp^ 289-316, 19F?. 



5. Mirasol. N. M. , "The Output of an M/G/ Queueing 
System is Poisson", Opns . Res . , 11, pp. 282-284, 1963. 



6. Lewis, P. A. W. and Uribe, L. , "The New Naval 
Postgraduate School Random Number Package LLRANDOMII , 
NPS 55-81-005 . Monterey, California, 1981. 



7. Quenouille, M. H. , "Notes an Bias in Estimation", 
Biometrika . 4, pp. 353-360, 1956. 



8. Tukey, J. W. , "Bias and Confidence in Not Quite Large 
Samples ', Ann . Math . Statist . , 29, p. 614, 1958. 



9. Miller, R. G. , "The Jackknife - a Review", Biometrika . 
61, pp. 1-15, 1974. 



10. Efron, B. , "Bootstrap Methods : Another Look at the 

Jackknife", Ann . Statist . , 7, pp. 1-26, 1979. 



11. Efron, B. , The Jackknife . the Bootstrap . and other 
Resampling Plans "" Technical Report No. 63, Dept' 6T 
Statistics, Stanford University, 1980. 



12. Feller, W. An Introduction to Probability Theory and 
its Applications " Vol. T] New~York, Wiley, 1966. 



13. Efron, B. , "Nonparametric Standard Errors and 
Confidence Intervals . Canadian Journal of Statistics. 
9, pp. 139-172, 1981. 



71 



INITIAL DISTRIBUTION LIST 



No. Copies 

1. Defense Technical Information Center 2 

Cameron Station 

Alexandria, Virginia 22304-6145 

2. Library, Code 0142 2 

Naval Postgraduate School 

Monterey, California 93943-5002 

3. Professor Patricia A. Jacobs, Code 55Jc 1 

Department of Operations Research 

Naval Postgraduate School 
Monterey, California 93943-5000 

4. Professor Donald A. Gaver, Code 55gv 2 

Department of Operations Research 

Naval Postgraduate School 
Monterey, California 93943-5000 

5. Professor Alan W. Washburn, Code 55 1 

Department of Operations Research 

Naval Postgraduate School 
Monterey, California 93943-5000 

6. CPT. Choi, Eui Joung 1 

Devision of System Analysis 

P. 0. Box 201-17 

Youngdeongpo Gu Singil 7 Dong 
Seoul Korea 

7. LT. Park, Dong Keun 10 

Incheon si, Booku Sangok Dong 307 

Bupyeong Hyeun Dai Apt 113-103 
Seoul Korea 

8. LCDR. Kim, Won Bong 1 

458 Ramona Ave. #1 

Monterey, California 93940 



72 




h-rv- I O "I 



Thesis 

P157404 

c.l 



Park 

Parametric and non- 
parametric estimation 
of the mean number of 
customers in service 
for AN M/G/00 queue . 



Thesis 
P157404 
c. 1 



21752G 



Park 

Parametric and non- 
parametric estimation 
of the mean number of 
customers in service 
for AN M/G/00 queue. 



